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Abstract: Soil organic carbon (SOC) is an effective indicator of soil 
fertility and productivity, and it varies spatially and temporally in relation 
to other soil properties. Spatial variability of SOC in the forestlands of 
northeast China was characterized using geostatistics. Soil samples at the 
depths of 0-20 cm, 20-40 cm and 40-60 cm were collected from six¬ 
ty-three temporary plots to evaluate SOC concentration and density 
(SOCD) and other soil properties. We analyzed correlations between 
SOC and soil properties. Soil organic carbon concentrations were high. 
The total amount of C stored in soil (0-60 cm) was 16.23 kg-m' with the 
highest SOCD of 7.98 kg-m' 2 in topsoil. Soil properties in most cases 
differed by horizon, suggesting different processes and effects in each 
horizon. Soil organic carbon had positive relationships with total N, P 
and K as well as readily available K, but did not show a significant posi¬ 
tive correlation with available P. Spatial factors including elevation, 
slope and aspect affected SOC distribution. Soil organic carbon at 0-60 
cm had strong spatial autocorrelation with nugget/sill ratio of 5.7%, and 
moderate structured dependence was found at 0-20 cm, which indicated 
the existence of a highly developed spatial structure. Spatial distributions 
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of SOC concentration and SOCD were estimated using regres- 
sion-kriging, with higher prediction accuracy than ordinary kriging. The 
fractal dimension of SOC indicated the preferential pattern of SOC dis¬ 
tribution, with the greatest spatial heterogeneity and strongest spatial 
dependence in the northeast-southwest direction. 

Keywords: northeast China, soil organic carbon, spatial variability, 
geostatistics 

Introduction 

Soil organic carbon (SOC) is a key determinant of soil quality 
(Stocking 2003) and represents one of the largest reservoirs of 
organic carbon in the global carbon cycle that can influence 
global warming (Lutzow et al. 2006). In terrestrial ecosystems, 
forest soil carbon is the largest carbon sink. Soil organic carbon 
storage in global forest ecosystems was estimated as 
402x10 12 -787x10 12 kg, or 25%-50% of the global soil carbon 
storage (Arrouays and Pelissier 1994). In northeast China, con¬ 
siderable organic carbon is stored in soils because of the dense 
forest vegetation and the thick accumulation of organic matter. 
Humic cambisols (Dark brown forest soils, in Chinese soil tax¬ 
onomy) are typical soils mainly found in warm-temperate, mixed 
coniferous and broad-leaved forests and cold-temperate conifer¬ 
ous forests in northeast China. The soils were higher in natural 
fertility and supported dense vegetation. The soils of 40.2x10 6 ha 
are mainly located in the provinces of Heilongjiang, Jilin and 
Inner Mongolia (Zhao et al. 2004). These areas have been, 
therefore, among the most important for timber production and 
as a sink for organic matter on a global scale. Humic cambisols 
are also considered to store substantial amounts of organic matter, 
thereby serving as C0 2 sinks. 

Geostatistics is an effective method for study of spatial distri¬ 
bution characteristics and their variability. Geostatistical meth¬ 
ods quantify spatial distributions based on the spatial scale of the 
study area, distance between detecting points, and spatial pat¬ 
terns of modeling semivariograms. They have been widely ap- 
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plied to evaluate spatial correlation in soils and to analyze the 
spatial variability of soil characteristics, such as soil physical and 
chemical properties (Vieira et al. 2007; Zheng et al. 2009) and 
biological properties (Fromm et al. 1993). Similar to other soil 
properties, SOC concentration and stock are highly variable over 
space (horizontally and vertically) and time (Wigginton et al. 
2000). Soil organic carbon density can vary due to intrinsic or 
extrinsic factors. Intrinsic variability is caused by natural varia¬ 
tion (Cambardella et al. 1994) and extrinsic variability results 
from external factors imposed on sites (Rao and Wagenet 1985). 
Spatial correlation of SOC can be explored by geostatistical 
analysis. On a national scale, Lacelle established the soil land¬ 
scape and SOC description database of Canada, and calculated 
that SOC storage at depths of 0-30 cm was 72.8xl0 12 kg 
(Lacelle 1998). Bernoux estimated 0-30 cm SOC pool of 
36.4(±3.4)xl0 12 kg in Brazil with soil and vegetation maps 
(Bernoux et al. 2002). Soil carbon reserves in different soil hori¬ 
zons were studied in New Zealand (Scott et al. 2002) and Congo 
(Schwartz and Namri 2002). In China, with the aid of GIS spatial 
analysis, soil organic carbon density (SOCD) and storage were 
estimated for the entire country and a high degree of spatial var¬ 
iability was documented. Many studies investigated the spatial 
variability of SOC on a regional scale. Xue et al. (2003) sur¬ 
veyed the Qaidam basin and calculated the SOCD to be 6.11 
kg-m" 2 . Using current 1:1000000 digital soil maps and 736 typi¬ 
cal soil profiles in northeast China, organic carbon storage in 
topsoil and its spatial distribution were estimated, and the aver¬ 
age SOCD at 0-1 m depths was 16.13 kg-m' 2 (Sun et al. 2004). In 
the tropical area of south China, SOC storage and spatial distri¬ 
bution were studied by combining soil data with a digital eleva¬ 
tion model (DEM). Similar studies were conducted in Jiangsu, 
Guangdong and Anhui provinces (Jiang et al. 2005; Gan et al. 
2003; Cheng and Xie 2009). According to previous studies, the 
initial application of GIS was to soil and land resource survey 
and evaluation, with focus on assessment of farmland soil deg¬ 
radation and mapping. Forestland, as a high potential land use 
type for carbon storage, has not yet gained much attention. 
Moreover, forest SOC storage or SOCD estimation on regional 
scales has large uncertainties because it is mainly based on soil 
survey data and literature and lacks detailed and reliable meas¬ 
urement. Ordinary kriging (OK) has been widely used for spatial 
interpolation in plains, but its predictive capability is limited in 
regions of complex topography with highly developed spatial 
conditions such as mountains and hills. Some ancillary variables 
such as terrain attributes and vegetation indices extracted from 
DEMs and remote sensing images are available for digital soil 
maps. Combining maps of soil types with soil properties infor¬ 
mation from point observations can improve spatial prediction of 
soil properties compared with prediction from a soil map only 
(Goovaerts and Journel 1995). 

The forest area in northeast China covers 3590xl0 6 ha, ac¬ 
counting for 24% of the total forest area of China, and the total 
forest stock volume amounts to 3.2x10 m (Department of For¬ 
est Resources Management 2010). It is the largest forest area and 
the transition region between broad-leaved mixed forest to the 
south and pure coniferous forest to the north. Research in this 
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region has mainly focused on the relationship between forest and 
climate change and on spatial heterogeneity of organic carbon 
storage in vegetation layer. However, the forest soil carbon pool 
has seldom been studied and most spatial interpolation used OK. 
Information on the spatial distribution of SOC is important for 
evaluating soil ecological functions and understanding soil car¬ 
bon sequestration processes. The use of geostatistics in combina¬ 
tion with DEM and regression-kriging (RK) can provide such 
information and help to define and improve the accuracy of rep¬ 
resentations of the distribution of soil parameters across land¬ 
scapes. Increased accuracy in the delimitation of SOC content 
classes could assist in meaningful soil quality prediction and 
assessment of the impacts of soil redistribution. The objectives 
of the present study were to evaluate the carbon-related soil 
properties and terrain attributes, and to determine the spatial 
variability and SOC distributions in different soil horizons of 
Humic cambisols in forests of northeast China using geostatisti¬ 
cal analysis. This information will provide a scientific basis for 
sustainable forest management and ecological environment im¬ 
provement in this region. 

Materials and methods 

Study area 

The study was carried out in Jincang Forest Farm, a typical forest 
area in the north temperate region of Wangqing County, Jilin 
Province, China (Fig. 1). The study area was in the eastern part 
of the low-middle hill region of the Changbai Mountain range 
(43°19-43°23' N, 130°26'-130°37' E) at elevations of 562 to 
950 m. The study area has a north temperate, monsoon climate 
with hot summers, cold winters, and abrupt changes of tempera¬ 
ture. Average annual air temperature is 3.9°C. Average annual 
precipitation is 547 mm, concentrated from July to September 
(1970-2010). Soils are characterized by thick and black topsoil, 
loamy texture, an acidic pH range, large amount of humus and 
high natural fertility. The predominant soils in the area are Hu¬ 
mic cambisols, together with meadow, peaty, swamp soil in 
some valleys that are developed from residual deposits of parent 
material such as granite, gneiss and basalt. Due to special land- 
forms, abundant precipitation and changeable climate, rich vege¬ 
tation resources are formed, including coniferous forest, 
broad-leaved forest, and mixed coniferous and broad-leaved 
forest. The major arbor species are Changbai larch ( Larix olgen- 
sis ), spruce ( Picea koraiensis ), Korean pine ( Pinus koraiensis 
Sieb. et Zucc.), maple ( Acer mono Maxim.), ash ( Fraxinus 
mandshurica Rupr .), Asian white birch ( Betula platyphylla Suk.) 
and linden ( Tilia amurensis Rupr.). 

Soil sampling and analysis 

In August 2011, a total of 63 temporary circular sample plots 
(250 m 2 each plot) were set up in the study area. Within each plot, 
a soil profile was dug to describe profile characteristics and sam¬ 
ples from the “bottom up” were taken with a ring sampler for 
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measurement of soil bulk density. The surface litter was removed 
prior to composite soil sampling. Soil samples at depths of 0-20, 
20-40, and 40-60 cm were collected with 5 soil cores each well 
mixed into a composite soil sample, which was further divided 
into two sets of sub-samples. One set of the sub-samples was air 
dried and passed through a 2 mm sieve for analysis of soil tex¬ 
ture, pH and available P and K; the remaining sub-sample was 
air-dried, homogenized and passed through a 0.25 mm sieve for 
analyses of SOC, total N, P and K (Li et al. 2011). Gravel and 
stones (fractions >2 mm) were sorted out and weighed. Soil tex¬ 
ture was measured using the hydrometer method (Wilde et al. 
1972). Soil pH was measured (1:2.5 = w/v) with a pH-meter 
(pH/ISE meter, model 710A; Orion). Soil organic carbon (SOC) 
concentration was determined by Walkley-Black wet oxidation 
method (Nelson and Sommers 1996; Bao 2000). Soil total N (TN) 
concentration was determined by Kjeldahl’s azotometer and total 
P (TP) by spectrophotometer, following wet digestion in concen¬ 
trated H 2 S0 4 (Bremner 1996). Total K (TK) was determined by 
flame spectrophotometer following wet digestion in HF-HC10 4 
(Jackson 1958; Knudsen et al. 1982). Available P (AP) was ex¬ 
tracted with HCI-NH 4 F and estimated by spectrophotometer 
(Bray and Kurtz 1945). Readily available K (AK) was extracted 
with 1.0 M NH 4 OAc and quantified by flame spectrophotometer 
(Carson 1980). Soil bulk density was determined by oven-drying 
the soil samples collected with a ring sampler at 105°C for more 
than 8 h. 



Fig. 1 Study area location and sampling design. 

In addition to soil sampling, the longitudes and latitudes of the 
sampling points were recorded with global positioning system 
(GPS). Forest type and survey, including environmental infor¬ 
mation, vegetation and management, for each sampling point 
were also recorded. A digital elevation model (DEM) of the 
study area was acquired at 5-m spatial resolution. The DEM and 
sampling points were geo-referenced to the Gauss Kruger coor¬ 
dinate system with Pulkovo 1942 GK Zone 22 (Fig. 1). 

Calculation of soil organic carbon density (SOCD) 

Soil organic carbon density (SOCD) is the organic carbon con¬ 
tent per unit area. We calculated SOCD in the profile with k 


layers ( S h kg-m' 2 ), on a ground-area basis up to 60 cm depth as 
follows: 

k k 

s k = 2A = Z c < x D < x E > x 0 - G ,)/i°o (!) 

1=1 i=1 

where, C z is SOC concentration (g-kg' 1 ), D { is soil bulk density 
(g-cm' 3 ), Ej is soil thickness (cm) and G, is volumetric percentage 
of the fraction >2 mm (rock fragments). S f is soil organic carbon 
density in the z'th horizon. In the present study, k= 3. 

Statistical analysis 
Descriptive statistics 

A descriptive statistical analysis was carried out to determine the 
mean, minimum, maximum, standard deviation (SD) and coeffi¬ 
cient of variation (CV) of the original data variables. The corre¬ 
lations among SOC, other soil properties and topographic factors 
such as elevation and aspect were also tested using correlation 
analysis and multiple linear stepwise regression. The significance 
of the difference between means of each property variable for 
different soil horizons was analyzed using one-way ANOVA. 
The statistical software SPSS 18.0 was used for statistical analy¬ 
sis. 

Geostatistical analysis 

The terrain analysis was based on a 5-m resolution digital eleva¬ 
tion model (DEM) generated from a 1: 5000 contour topog¬ 
raphical map. The elevation, slope in degrees, and aspect ex¬ 
pressed in positive degrees from 0 to 360 measured clockwise 
from north of the study site were generated from the DEM using 
ArcGIS 9.3 (Esri Company, RedLands, USA). Geostatistical 
analysis was performed using the software GS+ and ArcGIS 9.3. 
A sample semivariogram was calculated for each soil variable as 
follows (Goovaerts 1999; Wang and Shao 2013): 

y(h) = --- V N{h) [z(x j ) - z(x. Eh)] (^) 

2N(h)^ i=l 1 1 

where, z(x ; ) and z(x z +/z) are sample values at two points sepa¬ 
rated by the distance interval h , and N (h) is the total number of 
sample pairs for the lag interval h. For irregular sampling, it is 
rare for the distance between the sample pairs to be exactly equal 
to h. Therefore, h is often represented by a distance interval, y (h) 
is the sample semivariance for the distance lag h , which reflects 
spatial relationship between neighbor sampling points. 

The sample semivariogram was calculated and then fitted with 
a suitable theoretical model chosen from spherical, exponential, 
Gaussian, MateTn or Stein’s MateTn (Cressie 1993). The 
semivariogram models provide information about the spatial 
structure as well as the input parameters for kriging interpolation. 
Spatial dependence of soil variables can be evaluated by three 
important parameters: nugget (C 0 ), partial sill (Ci) and range in 
the semivariogram plot. The nugget is defined as the variability 
at a scale smaller than the sampling interval and/or sampling and 
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analytical error; the partial sill reflects the amount of spatial 
structural variance; and the range expresses the distance at which 
the semivariogram stabilizes around a limiting value. The spatial 
variation distance increasing to a constant value (sill) is known 
as range of spatial dependence (Cambardella et al. 1994), and it 
is affected by the sampling scale. The distribution distances of 
sampling points closer to the range are spatially correlated and 
the distances greater than the range are independent. In addition, 
the fractal dimension could be used to detennine the spatial 
variability of soil properties. When the fractal dimension values 
are small, soil properties have weak spatial dependence and sim¬ 
ple distribution, or vice versa (Burrough 1981, 1983; Wang and 
Shao 2013). 

Regression-kriging (RK) uses traditional regression methods 
that model the trend or drift in conjunction with geostatistical 
methods that model the random residuals from the regression 
(Hengl et al. 2004; Rivero et al. 2007). In RK the trend and 
residuals from the regression are fitted separately and then 
summed afterwards. It considers influence of geographic factors 
in the complex topography and simulates both spatial 
distribution trend and uncertainty. The trend item was separated 
by establishing regression equations between the auxiliary 
variable and target variable, and interpolation was applied to the 
residual with ordinary kriging (OK). The predicted value of the 
target variable was then obtained by spatial overlay between the 
trend item of regression prediction and OK value of the residual. 

Accuracy evaluation 

The prediction accuracy could be tested by comparing observed 
values with predicted values for soil nutrients, soil organic car¬ 
bon (SOC) and soil organic carbon density (SOCD) from sample 
plots. Mean prediction error ( M) and root mean square error ( R ) 
of prediction were adopted as the evaluation indices, and R was 
used to analyze the relative accuracy improvement value by 
comparing regression-kriging (RK) with ordinary kriging (OK) 
(Lai 2000; Kucharik et al. 2001). 



where, n is a sample point of test set, Z oi and Z pi are, respectively, 
observed values and predicted values, R 0 k is root mean square 
error of prediction using OK, and R rk is the prediction using RK. 
The prediction accuracy of RK is higher than OK when the R { 
value is positive, and the larger the value, the greater the accu¬ 
racy. In this study, 10 soil sampling points were randomly ex¬ 
tracted from 63 points to test the prediction accuracy of the re¬ 
maining 53 points. 
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Results 

Descriptive statistics of soil organic carbon (SOC) and other soil 
properties 

Table 1 shows the descriptive statistics of soil organic carbon 
(SOC) concentration and soil organic carbon density (SOCD). 
The mean of soil characteristics indicated the central tendency, 
and standard deviation (SD), coefficient of variation (CV) and 
skewness were used as the estimates of variability from each site. 
The average SOC concentration was 42.32 g-kg' 1 in the topsoil, 
and it decreased with depth until reaching 14.35 g-kg' 1 at 40-60 
cm. The CV was low for SOC concentration in 0-20 cm 
(CV<10%), indicating low variability of SOC in topsoil. How¬ 
ever, at soil depths of 20-40cm and 40-60 cm, SOC concentra¬ 
tion showed high CVs (67% and 80%). Accordingly, the SOCD 
was 7.98 kg-m' 2 for the topsoil, and it declined gradually with 
increasing depth to 5.23 kg-m' 2 at 20-40 cm and 3.06 kg-m' 2 at 
40-60 cm. Although there was a decreasing trend with increas¬ 
ing depth, the subsoil also contained a considerable amount of C, 
and the total amount of organic C stored to the depth of 60 cm 
was 16.23 kg-m' 2 , with a CV of 52%. 


Table 1: Descriptive statistics of soil organic carbon (n=6 3) 


Item 

Soil 

depth 

(cm) 

Mean 

SD 

Min. 

Max. 

CV 

(%) 


0-20 

42.32 

3.19 

11.63 

130.00 

7.54 

SOC concen- 

20-40 

24.49 

16.30 

8.71 

77.59 

66.56 

tration (g-kg’ 1 ) 

40-60 

14.35 

11.47 

2.97 

46.81 

79.93 


0-60 

27.40 

20.88 

2.97 

130.00 

76.20 


0-20 

7.98 

3.19 

1.35 

17.42 

39.97 

SOCD (kg-m’ 2 ) 

20-40 

5.23 

3.06 

0.86 

14.56 

58.51 

40-60 

3.02 

1.92 

0.42 

8.20 

63.58 


0-60 

16.23 

8.41 

7.63 

40.87 

51.82 


Notes: SOC is soil organic carbon; SOCD is soil organic carbon density; SD 
is standard deviation; CV is coefficient of variation. 

The descriptive statistics for other soil properties in the study 
area are shown in Table 2. Soil pH had a mean of 5.41-5.57 with 
similar SDs of 0.26-0.30 in three soil layers. The CVs of total N 
(TN), total P (TP) and available P (AP) all exceeded 50%, which 
indicated considerable spatial variation in these properties within 
the study area. The concentrations of TN, TP and AP had means 
of 1.07 g-kg' 1 , 0.38 g-kg' 1 and 25.49 mg-kg" 1 , respectively, at 
0-60 cm depth for all sampling points, and means and SDs var¬ 
ied by soil horizon. Mean values for all of these properties 
peaked at 0-20 cm and declined with increasing depth. 

Soil properties other than AP varied significantly by soil hori¬ 
zon. Concentrations of total N (TN), total P (TP), total K (TK) 
and readily available K (AK) declined with increasing depth in 
the following rank order: 0-20 cm > 20-40 cm > 40-60 cm with 
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significant differences between each horizon (p<0.05); this trend variation between horizons and better nutrient status in the top- 

was, however, not followed by available P (AP). In general, the soil, 

typical distribution of soil chemical properties suggested higher 


Table 2: Descriptive statistics of soil properties in different soil horizons 


Soil properties 

pH 

TN 

(g-kg 1 ) 

TP 

(g-kg 1 ) 

TK 

(g-kg 1 ) 

AP 

(mg-kg' 1 ) 

AK 

(mg-kg 1 ) 

Clay content 

(%) 

Gravel content 

(%) 

0-20 cm 










Mean 

5.41 a 

1.84 a 

0.51 a 

27.29 a 

22.27 

115.89 a 

28.94 a 

11.72 


SD 

0.26 

1.04 

0.25 

3.59 

11.48 

45.85 

2.13 

7.13 


Min. 

4.95 

0.57 

0.11 

17.27 

3.47 

80.00 

21.06 

0 


Max. 

6.22 

7.05 

1.10 

35.87 

46.69 

220.00 

33.44 

36.46 


CV (%) 

5 

57 

49 

13 

52 

40 

21 

54 


Skewness 

0.505 

2.998 

0.168 

-0.153 

0.107 

0.266 

0.233 

3.21 


20-40 cm 










Mean 

5.46 ab 

1.01 b 

0.40 b 

25.16 b 

19.79 

81.13 b 

25.93 a 

24.45 


SD 

0.27 

0.77 

0.23 

3.35 

11.54 

36.58 

1.78 

8.61 


Min. 

5.01 

0.21 

0.25 

16.81 

3.47 

30.00 

21.14 

3.31 


Max. 

6.09 

4.37 

1.18 

32.35 

41.99 

160.00 

31.67 

45.23 


CV (%) 

5 

76 

58 

13 

58 

45 

31 

60 


Skewness 

0.404 

2.258 

1.307 

-0.144 

0.163 

0.474 

-0.078 

1.570 


40-60 cm 










Mean 

5.57 b 

0.47 c 

0.29 c 

23.43 c 

18.94 

61.05 c 

21.56 ab 

27.73 


SD 

0.30 

0.31 

0.16 

3.47 

12.75 

27.30 

2.60 

10.12 


Min. 

4.89 

0.11 

0.09 

12.24 

4.13 

20.00 

17.79 

10.87 


Max. 

6.46 

1.75 

0.90 

30.73 

40.81 

140.00 

29.21 

60.55 


CV (%) 

5 

66 

55 

15 

67 

45 

29 

66 


Skewness 

0.601 

1.976 

1.198 

-0.848 

0.413 

0.819 

0.218 

3.61 


Notes: Different letters following the mean values indicate a significant difference (p < 

parison. TN, total N; TP, total P; TK, total K; AP, available P; AK, readily available K. 

0.05) between different horizons using 

a nonparametric multivariate com- 


Correlation analysis 

The correlation matrix of selected soil properties was shown in 
Fig. 2. Soil organic carbon (SOC) concentration was signifi¬ 
cantly and positively correlated with total N (TN), total P (TP), 
total K (TK) and readily available K (AK), but not with available 
P (AP). Soil TN was positively related with SOC concentration 
(r = 0.795, /><0.01), but its correlation with AP was weak. Cor¬ 
relations of AP with TN and TP were weak, while correlation of 
AP with AK became stronger with increasing depth. In each 
horizon, there were strong positive correlations between soil 
properties. In addition, TN, TP, TK and AK had significant posi¬ 
tive correlations with each other. The difference in SOC concen¬ 
tration was caused by natural organic matter distribution. To a 
certain extent, soil nitrogen levels affect SOC concentration 
through soil microorganisms. Different concentrations of P and 
K can be caused by different effects of organic matter, different 
pH values, and mineral composition. 

To the depth of 60 cm, the significant correlations between 
soil organic carbon (SOC) concentration and topographic attrib¬ 
utes indicated that SOC concentration tended to be greater at 
higher elevation. The correlation between SOC concentration 
and elevation was very significant and positive (r =0.532, 
7 ?<0.01). A very significant positive correlation was also found 
between soil organic carbon density (SOCD) and elevation (r = 


0.565, p<0. 01), but there was a very significant negative correla¬ 
tion between SOCD and slope (r = -0.451,/?<0.01). 



Fig. 2 : Correlation matrix between soil properties. **significant at prob¬ 
ability level of 0.01; ^significant at probability level of 0.05. TN is total N; 
TP is total P; TK is total K; AP is available P; AK is readily available K; 
SOC is soil organic carbon. 
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There were very significant correlations between SOC con¬ 
centration and SOCD (Table 3), and between these variables and 
cosine of slope aspect (p<0.01). SOC had, however, weak corre¬ 
lation with sine of slope aspect (which showed the degree of 
facing north), indicating that the more the aspect faced north the 
larger the SOC concentration. It tended to be greater on shady 
slopes at higher elevations. SOC in the topsoil had the same 
correlations with elevation and cosine of aspect, but this was not 
the case at depths of 20-40 cm and 40-60 cm. 


Table 3: Correlations between SOC concentration as well as SOCD and 
spatial factors («=63) 


Item 

Soil depth 

(cm) 

Elevation 

(m) 

Slope(°) 

cosa 

sina 


0-60 

0.532** 

-0.252 

0.409** 

-0.175 

SOC con- 


0-20 

0.514** 

0.178 

0.470** 

-0.052 

centration 


20-40 

0.173 

0.115 

0.055 

0.278 

g-kg% _ 


40-60 

0.184 

-0.318 

-0.212 

-0.010 


0-60 

0.565** 

-0.451** 

0.473** 

-0.018 

SOCD 

0-20 

0.452** 

-0.096 

0.375** 

-0.310 

(kg-m' 2 ) 

20-40 

0.193 

-0.105 

0.008 

0.244 


40-60 

0.219 

-0.365* 

0.204 

0.032 


Notes: ^significant at probability level of 0.01; ^significant at probability 

level of 0.05. SOC is soil organic carbon; SOCD is soil organic carbon den¬ 
sity. cosa=cosine of aspect; sina=sine of aspect. Aspect represents the degree 
of facing north. 

Multiple linear stepwise regression of soil organic carbon 
(SOC) concentration and soil organic carbon density (SOCD) on 
topographic factors yielded the following equations. Elevation 
and cosine of aspect (cosa) emerged as the optimal predictors of 
SOC and SOCD. 

Equations for 0-60 cm were: 


S OC =0.055 H+ 5.898 cosa-19.747 

7? 2 =0.344 (/?<().() 1) (6) 

V>cd=0.021 H+ 3.823 cosa-2.81, /? 2 =0.385 (p<0.01) (7) 

where, S oc is soil organic carbon concentration, S 0 C d is soil or¬ 
ganic carbon density, and H is elevation. The coefficients of 
determination for fitting SOC and SOCD were 0.344 and 0.385, 
respectively, y?<0.01 in both cases indicating good fits for both 
equations (6) and (7). 

Equations for 0-20 cm were: 

Sqc =0.119/7+13.787 cosa-52.806, 7? 2 =0.264 (p<0.01) (8) 

Vcd =0.020/7+2.224 cosa-7.143, 7? 2 =0.304 (p<0.01) (9) 

Geostatistical analysis 

Spatial variability of soil organic carbon (SOC) and semivario- 
gram models 

Geostatistical parameters of soil organic carbon (SOC) are 
shown in Table 4. After interpolating the residual with ordinary 
kriging (OK), the log-transformed data of SOC concentration at 
0-20 cm depth was fitted well by a Gaussian model having a 
nugget/sill ratio (N s r) of 5.7%, a range of 0.294 km, and very 
strong spatial autocorrelation. The large spatial range indicated a 
strongly structured regional pattern of SOC. The Gaussian model 
was best fitted to sample semivariograms of SOC at 0-60 cm 
depth with moderate structured dependence (N s r =47.8%). As 
for soil organic carbon density (SOCD), the A S r at 0-20 cm was 
15.4%, suggesting a considerable development of the spatial 
structure and high spatial correlativity. SOC at 0-60 cm depth 
was fitted by an exponential model and the spatial structure was 
moderate to well developed (A S r=49.3%). 


Table 4: Semivariogram models for SOC concentration and SOCD’s residual error 


Variable 

Data Item 

Model 

Nugget (C 0 ) 

Sill (Co+Ci) 

V SR (Co/Co+Ci) 

Range (km) 

SOC concentration (0-60 cm) 

residual error 

Gaussian 

190.70 

398.90 

0.478 

2.875 

SOC concentration (0-20 cm) 

residual error 

Gaussian 

4.90 

84.50 

0.057 

0.294 

SOCD (0-60 cm) 

residual error 

Exponential 

195.20 

412.40 

0.493 

3.170 

SOCD (0-20 cm) 

residual error 

Exponential 

2.84 

18.47 

0.154 

0.300 


Notes: SOC is soil organic carbon; SOCD is soil organic carbon density; Ns r is nugget/sill ratio; C 0 is Nugget; Ci is Sill. 


Accuracy evaluation for different prediction methods 
The regression-kriging (RK) method considers environmental 
factors and separates residual error to increase prediction accu¬ 
racy. Fifty-three sampling points were randomly extracted from 
63 points to conduct ordinary kriging (OK) interpolation for 
regression prediction residual error of soil organic carbon (SOC). 
Fig. 3 shows the spatial distribution of SOC concentration and 
soil organic carbon density (SOCD) interpolated using RK. At 
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0-60 cm, the spatial distribution of SOC mainly concentrated on 
northerly aspects. Taking the digital elevation model (DEM) in 
comparison, SOC concentration and SOCD were relatively rich 
at high elevations. Soil organic carbon at 0-20 cm showed that 
spatial distribution of SOC concentration had a number of 
high-value centers but was mainly concentrated in high elevation 
regions. The results suggest that slopes of northerly aspect served 
as sinks of organic matter and became the largest sources of 
SOC. 
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Fig. 3: Spatial distributions of soil organic carbon (SOC) and soil or¬ 
ganic carbon density (SOCD) interpolated by regression-kriging. 


The remaining 10 sampling points were used to test the accu¬ 
racy of the two prediction methods. Regression-kriging (RK) 
prediction was much more detailed on spatial variation and more 
closely approximated spatial distribution condition of the actual 


value (Table 5). Positive R\ values indicated higher accuracy 
using the RK method. 

Table 5: Comparison of accuracy using different prediction methods 


Variable 

OK 


RK 



M 

R 

M 

R 

R\r/o 

SOC (0-60 cm) 

1.13 

12.84 

0.68 

5.28 

58.88 

SOC (0-20 cm) 

1.01 

4.59 

-0.33 

3.91 

21.57 

SOCD (0-60 cm) 

0.32 

2.10 

0.17 

1.72 

18.09 

SOCD (0-20 cm) 

-0.96 

1.33 

-0.48 

1.02 

23.31 


Notes: M is mean prediction error, R is root mean square error of prediction, 
and R\ is the relative accuracy improvement value by comparing regres¬ 
sion-kriging (RK) with ordinary kriging (OK). 

The fractal dimension values of soil organic carbon (SOC) in 
different soil horizons are shown in Table 6. The greatest iso¬ 
tropic fractal dimension was found at 0-20 cm, indicating the 
most complex spatial pattern of SOC and the highest degree of 
uniformity. The fractal dimension at 20-40 cm was the smallest, 
indicating a relatively simple spatial pattern and stronger spatial 
dependence than other soil horizons. 

In terms of slope aspect, the greatest fractal dimension of soil 
organic carbon (SOC) at 0-20 cm was found in the north- 
west-southeast (NW-SE) direction, and the smallest in the 
northeast-southwest (NE-SW) direction, which had the greatest 
spatial heterogeneity and the strongest spatial dependence, and 
was also the direction of the longest transect through the study 
site (Fig. 3). As for the coefficient of determination, NE-SW 
direction showed a preferential pattern of SOC distribution. At 
20-40 cm, the NE-SW direction had the smallest fractal dimen¬ 
sion with the strongest spatial dependence; east-west (E-W) di¬ 
rection was the preferential pattern distribution for SOC. The 
fractal dimension of SOC at 40-60 cm depth was greatest in the 
NW-SE direction; and smallest in the NE-SW direction, which 
had greatest spatial heterogeneity. The coefficient of determina¬ 
tion showed that NE-SW had the preferential pattern of SOC 
distribution. 


Table 6: Fractal dimension of SOC in different soil layers 


Orientation 


0-20 cm 

20-40 cm 


40-60 cm 

Fractal Dimension 

Determination coefficients 
(R 2 ) 

Fractal 

dimension 

Determination 
coefficients (R 2 ) 

Fractal dimen¬ 
sion 

Determination coefficients 
(R 2 ) 

Isotropic 

1.874 

0.945 

1.781 

0.960 

1.812 

0.972 

S-N 

1.841 

0.827 

1.827 

0.815 

1.765 

0.933 

NE-SW 

1.810 

0.958 

1.723 

0.913 

1.737 

0.986 

E-W 

1.823 

0.756 

1.745 

0.940 

1.804 

0.652 

NW-SE 

1.854 

0.773 

1.778 

0.761 

1.836 

0.866 


Notes: SOC is soil organic carbon; S-N is south-north; NE-SW is northeast-southwest; E-W is east-west; NW-SE is northwest-southeast. 


g-kg' 1 and soil organic carbon density (SOCD) to the depth of 60 
Discussion cm was 16.23 kg-m' 2 , which is in agreement with previous stud¬ 

ies that estimated SOCD in a range of 16.55-27.02 kg-m' 2 in this 
The average soil organic carbon (SOC) concentration was 27.40 region (Zu et al. 2011; Sun et al. 2004). Because of differences in 
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the locations of study areas and management activities, SOC 
varies spatially. In addition, research methods are not standard¬ 
ized, and estimates are not free of uncertainty due to sampling 
depth and quantity. 

Soil organic carbon (SOC) and soil organic carbon density 
(SOCD) were higher in topsoil and declined with increasing 
depth. Distribution of SOC is related to vegetation types and 
succession, and human disturbance, so organic carbon in topsoil 
is greatly influenced by these factors. In the forest ecosystem, 
soil organic matter is mainly from litter and root decomposition. 
That surface soil can obtain more forest litter leads to the en¬ 
richment of SOC. The SOCD is associated with soil bulk density 
and stoniness (diameters >2 mm). Several studies found that 
SOCD increased with increasing soil depth (Yang and Wang 
2005). However, some studies in forests suggest that the changes 
in SOC stocks are due to the development of vegetation restora¬ 
tion and succession. For instance, vegetation restoration can 
increase non-active organic carbon concentration in surface soil 
(Su et al. 2005). Succession of restored vegetation can also in¬ 
crease SOC and nitrogen concentration (Christian et al. 2001). 

Correlation analysis quantifies the tendency of variables such 
as SOC, soil properties, and terrain attributes to vary in similar 
ways. Soil properties have indirect influence on SOC by 
affecting litter decomposition, microbial activity and organic 
carbon mineralization rates (Fayez and Esmael 2006). Soil total 
N (TN) concentration can promote the accumulation of biomass 
in understory vegetation and litter generation, thereby facilitating 
SOC accumulation. As a result, SOC, elevation, aspect, TN and 
total K (TK) were all interrelated, suggesting that the 
geographical characteristics and soil properties were correlated 
with SOC. Soil organic carbon (SOC) in all soil horizons had a 
positive relationship with TN, total P (TP), TK and readily 
available K (AK), in agreement with previous observations 
(Mladkova et al. 2005). A significant positive correlation 
between SOC and elevation indicated that SOC concentration 
tended to be greater at higher elevation. This is inconsistent with 
studies reporting greater SOC concentrations on the valley floor 
because of the deposition of water, solutes and sediments down 
along the slope (Mueller and Pierce 2003; Sumfleth and 
Duttmann 2008; Wang et al. 2001). Our findings agree with 
those of Tsuia et al. (2004) and Zhang et al. (2012) who 
documented increasing SOC with increasing elevation. 

Spatial patterns of soil organic carbon (SOC) suggest that or¬ 
ganic matter dynamics in the field were affected by topography. 
The spatial variability of SOC and soil organic carbon density 
(SOCD) differed by soil horizon. The spatial distribution of sur¬ 
face SOC demonstrated more developed spatial structure than the 
other two horizons, implying that topsoil might be more sensitive 
to environmental disturbances, whereas spatial variability of 
SOC at depths of 20-40 cm and 40-60 cm was greatly affected 
by random factors and suggested considerably weaker spatial 
correlativity. By combining the kriging map with the digital 
elevation model (DEM), the SOC concentration was related to 
elevation, indicating greater SOC in the higher terrain area. 
Similar findings with regard to the importance of elevation to 
SOC were reported by Geng et al (2011). However, some re- 
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search yielded the opposite results (Cheng et al. 2009; Zhang et 
al. 2012). This study showed that predictive capability of ordi¬ 
nary kriging (OK) was sometimes limited in regions of complex 
topography affected by intense cultivation and land use change. 
Recently, more studies have demonstrated regression-kriging 
(RK) can provide secondary information to improve the predic¬ 
tion accuracy (Lark and Webster 2006; Liu et al. 2010). It is 
possible that other useful ancillary information, such as soil and 
forest types, can be added in using RK to produce more accurate 
representations of the spatial distribution of SOC (Hengl et al. 
2004). The RK method with elevation as the predictor gave more 
details in space than the OK method in predicting spatial distri¬ 
butions of SOC variables. In the present study, it was difficult to 
use the spatial information reliably to establish correlations with 
the SOC variables by using RK analysis because study area 
spanned a small range of elevations. 

Conclusions 

Soil organic carbon (SOC) concentration and soil organic carbon 
density (SOCD) decreased significantly with increasing soil 
depth. Mean SOC concentration was 27.40 g-kg' 1 and SOCD in 
the whole soil profile (0-60 cm) was 16.23 kgnf . Correlations 
between SOC and soil chemical properties varied by soil depth. 
Soil organic carbon concentration showed significant positive 
relationships with total N (TN), total P (TP), total K (TK) and 
readily available K (AK), but not with available P (AP). A sig¬ 
nificant positive correlation was found between SOC and eleva¬ 
tion, as well as between SOC and slope aspect. Soil organic car¬ 
bon tended to be higher on sunny slopes. In the surface horizon, 
spatial variation produced by random factors was larger and had 
very small spatial correlation. With increasing depth, SOC was 
less sensitive to external influence with smaller spatial variation. 
The effect of terrain factors on SOC was very complex. Result¬ 
ing maps showed the spatial distribution of SOC in the study 
area based on the data of spatial dependency, and they can be 
used to delineate the areas with different impact of elevation. 
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